{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 05: weak imposition of Dirichlet BCs by a Lagrange multiplier (inf-sup computation)\n",
    "\n",
    "In this tutorial we compute the inf-sup constant of the saddle point problem resulting from the discretization of the following Laplace problem\n",
    "$$\\begin{cases}\n",
    "-\\Delta u = f, & \\text{in } \\Omega,\\\\\n",
    " u   = 0, & \\text{on } \\Gamma = \\partial\\Omega,\n",
    "\\end{cases}$$\n",
    "\n",
    "where $\\Omega$ is the unit ball in 2D, and for which the non-homogeneous Dirichlet boundary conditions are imposed by a Lagrange multiplier.\n",
    "\n",
    "The resulting eigenvalue problem is\n",
    "$$\n",
    "\\text{find } \\eta, u, \\lambda \\in \\mathbb{R} \\times V \\times M \\text{ s.t. }\\\\\n",
    "\\begin{cases}\n",
    "\\int_\\Omega \\nabla u \\cdot \\nabla v + \\int_\\Gamma \\lambda v = 0, & \\forall v \\in V,\\\\\n",
    "\\int_\\Gamma u \\mu = \\eta \\int_\\Gamma \\lambda \\mu, & \\forall \\mu \\in M\n",
    "\\end{cases}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "V = H^1(\\Omega),\\\\\n",
    "M = L^{2}(\\Gamma).\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mpi4py import MPI\n",
    "from petsc4py import PETSc\n",
    "from slepc4py import SLEPc\n",
    "from ufl import dx, grad, inner, Measure, TestFunction, TrialFunction\n",
    "from dolfinx import Constant, FunctionSpace\n",
    "from dolfinx.fem import locate_dofs_topological\n",
    "from dolfinx.io import XDMFFile\n",
    "from multiphenicsx.fem import assemble_matrix_block, DofMapRestriction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with XDMFFile(MPI.COMM_WORLD, \"data/circle.xdmf\", \"r\") as infile:\n",
    "    mesh = infile.read_mesh()\n",
    "    mesh.topology.create_connectivity_all()\n",
    "    boundaries = infile.read_meshtags(mesh, name=\"boundaries\")\n",
    "facets_Gamma = boundaries.indices[boundaries.values == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define associated measures\n",
    "ds = Measure(\"ds\")(subdomain_data=boundaries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenvalue problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function space\n",
    "V = FunctionSpace(mesh, (\"Lagrange\", 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define restrictions.\n",
    "dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)\n",
    "dofs_V_Gamma = locate_dofs_topological(V, boundaries.dim, facets_Gamma)\n",
    "restriction_V = DofMapRestriction(V.dofmap, dofs_V)\n",
    "restriction_V_Gamma = DofMapRestriction(V.dofmap, dofs_V_Gamma)\n",
    "restriction = [restriction_V, restriction_V_Gamma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define trial and test functions\n",
    "(u, l) = (TrialFunction(V), TrialFunction(V))\n",
    "(v, m) = (TestFunction(V), TestFunction(V))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define problem block forms\n",
    "a = [[inner(grad(u), grad(v)) * dx, - l * v * ds],\n",
    "     [- u * m * ds, None]]\n",
    "b = [[None, None],\n",
    "     [None, - l * m * ds]]\n",
    "b[0][0] = Constant(mesh, 0.) * inner(u, v) * dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assemble lhs and rhs matrices\n",
    "A = assemble_matrix_block(a, bcs=[], restriction=(restriction, restriction))\n",
    "A.assemble()\n",
    "B = assemble_matrix_block(b, bcs=[], restriction=(restriction, restriction))\n",
    "B.assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve\n",
    "eps = SLEPc.EPS().create()\n",
    "eps.setOperators(A, B)\n",
    "eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)\n",
    "eps.setDimensions(1, PETSc.DECIDE, PETSc.DECIDE)\n",
    "eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)\n",
    "eps.setTarget(1.e-5)\n",
    "eps.getST().setType(SLEPc.ST.Type.SINVERT)\n",
    "eps.getST().getKSP().setType(\"preonly\")\n",
    "eps.getST().getKSP().getPC().setType(\"lu\")\n",
    "eps.getST().getKSP().getPC().setFactorSolverType(\"mumps\")\n",
    "eps.solve()\n",
    "assert eps.getConverged() >= 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract leading eigenvalue\n",
    "eigv = eps.getEigenvalue(0)\n",
    "r, i = eigv.real, eigv.imag\n",
    "assert abs(i) < 1.e-10\n",
    "assert r > 0., \"r = \" + str(r) + \" is not positive\"\n",
    "print(\"Inf-sup constant: \", np.sqrt(r))\n",
    "assert np.isclose(np.sqrt(r), 0.125496)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
